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We present a modified bond- valence model of PbTiOs based on the principles of bond- valence and 
bond- valence vector conservation. The relationship between bond- valence energy and bond-order 
potential is derived analytically in the framework of a tight-binding model. A new energy term, 
bond-valence vector energy, is introduced into the atomistic model and the potential parameters 
are re-optimized. The new model potential can be applied to both canonical ensemble (NVT) and 
isobaric-isothermal ensemble (NPT) molecular dynamics (MD) simulations. This model reproduces 
the experimental phase transition temperature in NVT MD simulations and also exhibits the exper- 
imental sequence of temperature-driven and pressure-driven phase transitions in NPT simulations. 
We expect that this improved bond-valence model can be applied to a broad range of inorganic 
materials. 



The use of ferroelectric perovskite oxides in a vari- 
ety of technological applications has prompted extensive 
investigations of their structure and dynamics i 2 - First- 
principles density functional theory (DFT) calculations 
have played an important role in enhancing microscopic 
understanding of the relationships between composition, 
structure and properties £ — Despite the success of first- 
principles methods, the great computational expense 
and restriction to studying zero-temperature properties 
have driven the development of more efficient atom- 
istic and effective Hamiltonian potentials suitable for 
large-scale molecular dynamics (MD) simulations^— 

In particular, an atomistic potential based on the 
widely used bond- valence (BV) theory^ was developed 
by Shin et al*& BV-based atomistic potentials have since 
been used to study phase transitions^ and domain wall 
motion in PbTiOs 15 , as well as structure and dynam- 
ics in the classic 0.75PbMg 1/3 Nb 2/3 O3-0.25PbTiO3 re- 
laxor ferroelectric material^ The bond-valence theory, 
or bond- valence conservation principle, states that in a 
crystal structure, each atom i prefers to obtain a cer- 
tain atomic valence, Vb,i- The actual atomic valence 
Vi, for atom i can be obtained by summing over the 

bond valences Va, which can be calculated from an em- where <K r u) is a pair-wise repulsive potential depend- 



the general applicability of this type of atomistic po- 
tential. In addition, the potentials obtained in previous 
wor k 12 i 14 were found to be accurate for NVT simula- 
tions only, with incorrect ground state structures ob- 
tained when the constant volume constraint is lifted. In 
this paper, we show how the bond- valence potential can 
be derived from the second-moment of the local density 
of states (LDOS), extend the model to represent higher 
moments and show that this allows accurate simulations 
for both constant- volume and constant-pressure condi- 
tions. 

An analysis of the physics that gives rises to the bond- 
valence conservation principle shows that the bond- 
valence energy can be naturally derived from the second- 
moment bond-order potential, such as the well-known 
Finnis-Sinclair potential ] 19 i 2Q Within the framework of a 
tight-binding model 2 ^, the Finnis-Sinclair potential can 
be partitioned into atomic contributions as: 
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pirical inverse power-law relationship 17 ! 18 between bond 
valence and bond length r^: 
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7*0, and Cij are Brown's empirical parameters. The en- 
ergy contribution of the bond- valence is chosen to have 
the following form: 



ing on the distance between atom i and its nearest- 
neighboring atom j. The second term represents the 

(2) 

bonding energy; 7^ is a constant and ji\ ' is the second- 

(2) 

moment of the LDOS. The second moment ji\ mea- 
sures the width of the density of states distribution, 
and as shown by Cyrot-Lackmann and Ducastelle 2 ^— 
can be evaluated from summation over all the nearest- 
neighbor hopping paths that start and end on atom i: 
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where Si is a scaling parameter. 

Despite the success of the rather simple ten-parameter 
BV model potential! 2 - for PbTi0 3 , no rigorous quan- 
tum mechanical justification has been provided for the 
bond- valence potential energy, raising questions about 



where fyj is the averaged hopping integral between atom 
i and j. Because the overlap of atomic orbitals decays 



as exp(— kijrij) , Eq.(3) can be written as 
U FS 
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where n\ is the number of electrons at atom i, 
constant that scales the strength of the repulsive inter- 
actions between atom i and atom j, and bij scales the 
bonding interaction. The total energy will be minimized 
when the derivative of Ufs with respect to exp(— k^r^) 
is zero. This leads to 
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After rearranging the terms, we get the minimization 
condition to be 
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In the case of a two component system, there is only one 
type of atomic pair^ and therefore bij/aij is the same 
for all the atomic pairs. We can then rigorously define 
the interatomic bond- valence Vij as 
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where C = bfj/afj. It is evident that the energy mini- 
mization condition becomes 



(9) 



which is identical to the bond- valence summation rule. 
The bond valence of the individual bond Vij corresponds 
to the hopping integral with Vij oc /3fj . In the case of a 
more complicated ABO3 oxide, we would expect that a 
universal ratio of bij to a^- for all the cation-anion pairs 
is still a reasonable assumption. 

After substituting Eq.(8) into Eq.(5) and expanding 
in a Taylor series around Vb,i, we obtain: 
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Thus, the Finnis- Sinclair energy term (Eq.(3)) and the 
bond- valence energy (Eq.(2)) are identical for small de- 
viations away from Vo,i. 

Compared to the bond-order potential, the applica- 
tion of the bond-valence model does not require extra 
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FIG. 1: (Color online) Schematic representation of 
bond-valence vector summation around Ti in (a) cubic 
PbTiOs and (b) tetragonal PbTiOs. Gray, blue and 
red balls denote Pb, Ti and O. The back arrow scales 
the magnitude of the bond-valence. 



efforts to parametrize hopping integrals, because the 
bond-valence parameters for a wide variety of atomic 
paris are already known from crystallography. 13 Since 
the bond- valence model is a second-moment bond-order 
potential, its limitations, such as the inability to obtain 
the correct ground state structure in NPT simulations, 
are likely due to the fact that the second moment only 
accounts for the width of LDOS but does not reflect 
its shape. One consequence of this is that the BV 
energy depends only on the total valence and is entirely 
insensitive to the number of bonds or their relative 
strengths. This feature of second-moment models 
makes it difficult to distinguish between competing 
crystal structures, which are controlled by the higher 
moments. 20 Therefore, a systematic way to improve the 
bond-valence model is to include the contributions of 
higher moments of the LDOS (such as fourth- moment) 
to the total energy i 26 i 27 

We choose the bond- valence vector sum (BVVS) 13 i 28 
to reflect the change of the fourth moment of the LDOS. 
The bond-valence vector is defined as a vector lying 
along the bond with magnitude equal to the bond- 
valence, as shown in Figure 1. The changes in the local 
symmetry of the bonding environment that are reflected 
by the BVVS affect the LDOS shape and the value of the 
fourth LDOS moment. A simple argument is presented 
in the Appendix to illustrate the relationship between 
the fourth moment of the LDOS and the sum of the 
bond- valence vectors. 

For many materials, it has been shown that the 
ground-state structure favors symmetric local bonding 
environment and a zero BVVS. Therefore, the crite- 
rion of BVVS = for the ground-state structure has 
been suggested as a complement to the original bond- 
valence conservation principle . 13 ! 28 However, this is not 
true for crystal structures in which symmetry break- 
ing (BVVS 7^ 0) becomes significant due to electronic- 
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structure driven distortions, such as the second order 
Jahn- Teller distortion exhibited by Ti atoms in an oc- 
tahedral environment and the stereochemical long-pair 
driven distortion present in Pb 2+ and Bi 3+ cations. The 
BVVS can thus be considered as a measure of local sym- 
metry breaking. We therefore generalize this principle 
by proposing that each ion has a desired length of bond- 
valence vector summation. The bond- valence vector en- 
ergy, Ebvv, is defined as 



TABLE I: Optimized potential parameters of modified 
bond-valence model. Then angle potential parameter 
k is 3.14 meV/(deg) 2 









qp{e) Sp(eV) Dp 


Pb Ti O 


Vb }j 8 V ,/3 


Pb 1.38177 0.31646 2.23180 
Ti 0.99997 0.11888 
-0.79391 1.52613 


- 2.06058 1.70871 
1.29532 
1.88109 


2.00 0.40297 
4.00 0.46541 
2.00 



£wv = £A(Vf-V 2 / (11) 

i 

where Di is the scaling factor, is the calculated bond- 
valence vector summation and Vo,i is the desired value 
of bond-valence vector summation. The value of Vo,i 
can be computed using the optimized atomic positions 
in the lowest-energy structure from first-principles. 

The interatomic potential for our modified bond- 
valence model is given by: 

E = E c + E r + Ebv + Ebvv + E a (12) 

£ c = Ef^ (13) 

Kj J 

N oct 

^ = fc£(C+C+C) ( 15 ) 



where E c is the Coulomb energy, E r is the short-range 
repulsive Lennard- Jones energy. It is noted that all the 
atomic orbitals are approximated as s-type and only av- 
eraged hopping integrals between neighboring atoms are 
used in both Finnis- Sinclair potential and the BV en- 
ergy. However, bonding in PbTiOa involves p and d or- 
bital overlaps, which do display angular dependence. To 
introduce the dependence of energy on the interatomic 
angles, we used an angle potential, E a , which prevents 
unphysically large tilting of oxygen octahedra. The po- 
tential parameters required from fitting can be summa- 
rized as follows: spring constant fc, charges ^, scaling 
factors Si and Di for each species and short-range repul- 
sion parameters, Bij, for all the pairs of Pb-Ti, Pb-O, 
Ti-0 and O-O. Due to the charge- neutrality constraint, 
there are thirteen independent parameters. 

Figure 2 shows our parameterization protocol. The 
optimization of the potential parameters is performed 
using simulated annealing (SA) global optimization 
method to fit a database of structural energy differences 
and atomic forces (E & F) derived from ab initio DFT 
calculations with the Abinit code^ We used the 2x2x2 
supercell as the reference structure. The energy and 
atomic forces are computed with 2x2x2 Monkhorst- 
Pack fc-point mesh 30 using PBEsol 31 as the exchange- 
correlation energy functional. We start with an ini- 




FIG. 2: Potential optimization protocol used in this 
work. 



(a) 





0.5 




0.4 


em 


0.3 








0.2 


isp 


0.1 


Q 






1 — 1 — i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — r 



Pb 



<»-ooooooooooooo 



-ooooooooooooooooooo^-' 



J i I i I i I i I i I i L 



100 200 



300 400 500 600 
Temperature (K) 



700 800 900 




100 200 



300 400 500 600 
Temperature (K) 



700 800 900 



FIG. 3: (Color online) Temperature dependence of 
(a) atomic displacements of Pb and Ti, and (b) spon- 
taneous polarization obtained from NVT MD sim- 
ulations with lattice constants fixed to experimental 
values. 
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FIG. 4: (Color online) Temperature-dependent properties of PbTiOs obtained from NPT simulations, (a)-(e) 
Profiles of lattice constants from MD simulation, (f) Atomic displacements of Pb and Ti as a function of temper- 
ature, (g) Spontaneous polarization as a function of temperature. 



tial database that contains the lowest-energy tetrago- 
nal structure, strained tetragonal structures, the lowest- 
energy cubic structure, strained cubic structures, and 
randomly picked orthorhombic structures with various 
lattice constants. After each SA run, the optimized po- 
tential parameters are used to perform constant-stress 
MD simulations to generate equilibrium structures at 
various temperatures, which are then put back to the 
database. The process is continued until the energies 
of the structures sampled during MD simulations accu- 
rately reproduced and the simulation therefore repro- 
duces the DFT ground state and correct phase transi- 
tion behavior. 

Table I presents the optimized potential parameters. 
To account for the overestimation of the PbTiOs c/a 
ratio by PBEsol (c/a= 1.10 versus c/a= 1.07 experimen- 
tally) 32 , we adjusted Brown's empirical parameter rj- 
to make the Vp for Pb, Ti and O reach their atomic 
valences in the lowest-energy tetragonal structure. The 
value of preferred vector summation is then calculated 
with the modified empirical parameters. We found that 
the oxygen atoms do not have a preference for a spe- 
cific value of bond- valence vector sum. This is because 
in perovskites, some oxygen atoms are highly displaced 
(I Vo | > 0), while others stay around the high-symmetry 
point (| Vo| = 0). Therefore, we set the scaling factor 
Do fc> r oxygen atoms to be zero. During the parame- 
ter fitting, we also consistently found that the optimized 
scaling factor Sti is zero, which implies that for PbTiOs 
only two scaling factors of the bond- valence energy are 
independent, due to the bond- valence conservation prin- 



ciple. 

Using this optimized model potential, we studied the 
temperature dependence of lattice constants, polariza- 
tion and displacements of Pb and Ti using an 8 x 8 x 8 
supercell. We first performed canonical-ensemble MD 
simulation with lattice constants fixed to experimen- 
tal values, using the Nose-Hoover thermostat to con- 
trol the temperatures. For these simulations, we ob- 
tained 760 K for the ferroelectric-to-paraelectric first- 
order phase transition temperature T c , shown in Fig- 
ure 3. This agrees well with the experimental T c of 
765 K and is an improvement relative to the 550 K val- 
ues obtained by the NVT calculations with the origi- 
nal BV potential ] 12 i 14 We then used the new potential 
in NPT simulations, with the pressure maintained at 
0.1 MPa by the Parrinello- Rahman barostat. For the 
ground state structure at K, we obtained the lattice 
constant a=3.865 A and c/a= 1.12. The c/a ratio in MD 
is slightly larger than the PBEsol DFT value. Figure 
4 displays the temperature dependence of lattice con- 
stants, spontaneous polarization and atomic displace- 
ments of Pb and Ti obtained from NPT simulations. 
As temperature increases, the c/a ratio decreases grad- 
ually, together with the polarization and atomic dis- 
placements. The phase transition from tetragonal to cu- 
bic occurs at 400 K, lower than the experimental value. 
The rather large magnitude of spontaneous polarization 
and the atomic displacements at temperatures just be- 
low T c are due to the overestimated tetragonality of the 
PBEsol functional and the resulting potential. 
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ergy (£?dw) is calculated by 



(b) 




AW* + (c/a) 2 



N lC /y/l + {c/a) 2 



FIG. 5: (Color online) Simulated domain wall using 
modified bond- valence model, (a) 180° domain wall 
constructed with a 12x4x4 supercell; (b) 90° domain 
wall with Ni = 12, N 2 = 4, N 3 = 4. 
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FIG. 6: (Color online) Pressure-induced phase tran- 
sitions of PbTiOs obtained from MD simulations. 
Lattice axes coincide with the Cartesian axes (a along 
b along y and c along z). 



We find that the new potential is capable of describ- 
ing the domain wall (DW) energetics and structure. The 
supercell used to model the domain wall is constructed 
following the method in Ref. 32. The domain wall en- 



En — ^bu 



ik 



(16) 



where En is the energy of the supercell, £?buik is the 
energy of a single-domain supercell of the same size, 
and Sbw is the area of the domain wall. Figure 5(a) 
presents simulation of 180° Pb-centered domain walls at 
10 K. The computed domain wall energy is 161 mJ/m 2 
and agrees very well with 170 mJ/m 2 obtained via 
PBEsol DFT calculations (with a 8 x 1 x 1 supercell). 
To simulate a 90° domain wall, we used a supercell 
with N\ = 12, N2 = 4, and N$ = 4, as shown in 
Figure 5(b). The dimensions of the supercell are fixed 
to the values calculated based on experimental lattice 
constants of tetragonal PbTiOs. The domain wall 
energy is estimated to be 70 mJ/m 2 and also shows a 
satisfying agreement with the PBEsol DFT value of 64 
mJ/m 2 (with a 8 x 1 x 1 supercell). We note that the 
BV potential is highly efficient, as all the interactions 
are pair- wise and only depend on distance. This allows 
simulations of 30x30x4 supercell (18000 atoms) for 15 
ps with a 0.3 fs timestep using only 33678 seconds of 
CPU time (10 minutes of wall clocktime parallelized 
over 48 CPUs on the Cray XT5 supercomputer at the 
Navy DoD Supercomputing Resource Center). 

We have also examined the performance of the poten- 
tial in simulations of pressure-induced phase transitions 
in PbTi03 with a 10 x 10 x 10 supercell. Figure 6 
shows the pressure dependence of lattice constants and 
polarization. We found two phase transitions, at 7 GPa 
and 15 GPa. Below 7 GPa, the structure is ferroelectric. 
The tetragonality decreases with increased pressure 
and the magnitude of polarization along the long 
axis reduces accordingly. Above 7 GPa, the c/a ratio 
becomes 1 but the structure maintains ferroelectricity 
up to 15 GPa. The polarization is mostly along [111] 
direction between 7 GPa and 15 GPa, suggesting the 
existence of rhombohedral phase. The polarization 
disappears when the pressure exceeds 15 GPa and the 
structure becomes the centrosymmetric paraelectric. 
Our simulated results are consistent with Wu and 
Cohen's first-principles studie a 34 i 35 and recent exper- 
imental results by Ahart et alS^ We did not find any 
reentrance of ferroelectricity up to 60 GPa. 



We improved our bond- valence model by introducing 
a new energy term, bond-valence vector energy, based 
on an extension of the bond- valence vector conservation 
principle. This new 12-parameter potential reproduces 
the polarization, ferroelectric instability and phase 
transition temperature in NVT simulations, and also 
captures the temperature-driven phase transition quali- 
tatively in NPT simulations. Both calculated 180° DW 
energy and 90° DW energy using this new potential 
are in good agreement with DFT values. This new 
potential is efficient enough to simulate large domain 
walls. The studies of pressure-induced phase transition 
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with the new potential show two phase transitions, 
consistent with previous experimental studies. We hope 
that this improved bond-valence model will also be 
applied to other oxides due to its simplicity, efficiency 
and accuracy. 

S.L. was supported by the NSF through Grant 
CBET-0932786. H.T. was supported by the US DOE 
BES under Grant No. DE-FG02-07ER46431. LG. was 
supported by the Energy Commercialization Institute. 
T.Q. and A.M.R. were supported by the US ONR 
under Grant No. N00014-11-1-0578. Computational 
support was provided by the Center for Piezoelectrics 
by Design, and by the DoD HPCMO. 



APPENDIX 

The bond valence of an individual bond Vij has been 
shown to be proportional to the square of hopping in- 
tegral fyj. Both the bond- valence vector sum, V^, and 

the fourth-moment of the LDOS, /4 4 \ can re A ec t the 
change of local symmetry of bonding environment. Fig- 
ure Al gives an example of a one-dimensional AB alloy. 
The desired bond valence of A — B in the undistorted 
structure is set to be a, and therefore the hopping in- 
tegral is equal to ^/ja, where 7 is a constant. It is 
easy to calculate that the bond valence summation and 
/i^ 2 ) at atom A are 2a and 27a, respectively. The small 
displacement in the distorted lattice changes the bond- 
valence of the longer A — B bond to (a — 6) and the 
shorter one to (a + 5). Accordingly, the hopping inte- 
gral for the longer A — B is y/"y(a — 6) and the shorter 
one \Ay(a + 5). The bond-valence conservation princi- 
ple holds in both structures. However, the changes 
from zero in the undistorted structure to 25 in the dis- 
torted structure, and the is reduced from 6j 2 a 2 
to 67 2 a 2 — 2^ 2 S 2 . It is evident that only the hopping 
path involving the next-nearest neighbors contributes 
to the change of fourth- moment. Therefore, the change 
of fourth moment, A/^ 4 \ can be approximated with 
(|Vi| - IV^oI) 2 . We choose V? instead of |V*| in the 
formula of Ebvv to make sure Ebvv is a differentiate 
function for each VV 
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FIG. Al: Hopping integrals in one-dimensional AB 
alloy. Empty and filled circles represent elements A 
and B. 
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